\defmodule {InverseGaussianProcess}

The inverse Gaussian process is a non-decreasing process
where the increments are additive and are given by the
inverse gaussian distribution,
\externalclass{umontreal.iro.lecuyer.probdist}{InverseGaussianDist}.
With parameters $\delta$ and $\gamma$, the
time increments are given by
\externalclass{umontreal.iro.lecuyer.probdist}{InverseGaussianDist}$(\delta  dt/\gamma,
\delta^2 dt^2)$.


[We here use the inverse gaussian distribution
parametrized with {IGDist}$(\mu,\lambda)$, where $\mu=\delta/\gamma$
and $\lambda=\delta^2$.  If we instead used the parametrization
$\textrm{IGDist}^\star(\delta, \gamma)$,
then the increment distribution of our process would have been written
more simply as $\textrm{IGDist}^\star(\delta  dt, \gamma)$.]

The increments are generated by using
the inversion of the cumulative distribution function.
It therefore uses only one \externalclass{umontreal.iro.lecuyer.rng}{RandomStream}.
Subclasses of this class use different generating methods and some need
two \externalclass{umontreal.iro.lecuyer.rng}{RandomStream}'s.

The initial value of this process is the initial observation time.


\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        InverseGaussianProcess
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.stochprocess;\begin{hide}
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;

\end{hide}

public class InverseGaussianProcess extends StochasticProcess \begin{hide} {

    protected RandomStream       stream;

    protected double   delta;
    protected double   gamma;

    protected double   deltaOverGamma;
    protected double   deltaSquare;
    // mu and lambda are the common names of the params for InverseGaussianGen.
    protected double[] imu;
    protected double[] ilam;

    // Number of random streams needed by the current class
    // to generate an IG process.  For this class, = 1, for subclasses
    // will sometimes be, = 2.
    int numberOfRandomStreams;

\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}\begin{hide}
    // needed by InverseGaussianProcessMSH
   protected InverseGaussianProcess()  { }
\end{hide}

   public InverseGaussianProcess (double s0, double delta, double gamma,
                                  RandomStream stream) \begin{hide} {
        this.x0 = s0;
        setParams(delta, gamma);
        this.stream = stream;
        numberOfRandomStreams = 1;
    }\end{hide}
\end{code}
\begin{tabb} Constructs a new \texttt{InverseGaussianProcess}.
The initial value $s0$ will be overridden by $t[0]$ when
the observation times are set.
\end{tabb}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}\begin{hide}

   public double[] generatePath() {
        double s = x0;
        for (int i = 0; i < d; i++) {
            s += InverseGaussianGen.nextDouble(stream, imu[i], ilam[i]);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }\end{hide}

   public double[] generatePath (double[] uniforms01) \begin{hide} {
        double s = x0;
        for (int i = 0; i < d; i++) {
            s += InverseGaussianDist.inverseF(imu[i], ilam[i], uniforms01[i]);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }\end{hide}
\end{code}
\begin{tabb}
Instead of using the internal stream to generate the path,
uses an array of uniforms $U[0,1)$.  The array should be
of the length of the number of periods in the observation
times. This method is useful for \class{NormalInverseGaussianProcess}.
\end{tabb}
\begin{code}

   public double[] generatePath (double[] uniforms01, double[] uniforms01b) \begin{hide} {
       throw new UnsupportedOperationException("Use generatePath with 1 stream.");
    }\end{hide}
\end{code}
\begin{tabb}
This method does not work for this class, but will
be useful for the subclasses that require two streams.
\end{tabb}
\begin{code}\begin{hide}

   public double nextObservation () {
        double s = path[observationIndex];
        s += InverseGaussianGen.nextDouble(stream,
                  imu[observationIndex], ilam[observationIndex]);
        observationIndex++;
        observationCounter = observationIndex;
        path[observationIndex] = s;
        return s;
    }\end{hide}

   public void setParams (double delta, double gamma) \begin{hide} {
        this.delta = delta;
        this.gamma = gamma;
        deltaOverGamma = delta/gamma;
        deltaSquare    = delta*delta;
        init();
    }\end{hide}
\end{code}
\begin{tabb} Sets the parameters.
\end{tabb}
\begin{code}

   public double getDelta() \begin{hide} {
        return delta;
    }\end{hide}
\end{code}
\begin{tabb} Returns $\delta$.
\end{tabb}
\begin{code}

   public double getGamma() \begin{hide} {
        return gamma;
    }\end{hide}
\end{code}
\begin{tabb} Returns $\gamma$.
\end{tabb}
\begin{code}

   public double getAnalyticAverage (double time) \begin{hide} {
        return delta*time/gamma;
    }\end{hide}
\end{code}
\begin{tabb} Returns the analytic average which is $\delta t/ \gamma$,
   with $t=$ \texttt{time}.
\end{tabb}
\begin{code}

   public double getAnalyticVariance (double time) \begin{hide} {
        return delta*delta*time*time;
    }\end{hide}
\end{code}
\begin{tabb} Returns the analytic variance which is $(\delta  t)^2$,
   with $t=$ \texttt{time}.
\end{tabb}
\begin{code}\begin{hide}

    protected void init () {
        super.init(); // set path[0] to s0
        if(observationTimesSet){
            x0 = t[0];
            path[0] = x0;
            imu  = new double[d];
            ilam = new double[d];
            for (int j = 0; j < d; j++)
            {
                double temp = delta * (t[j+1] - t[j]);
                imu[j]  = temp / gamma;
                ilam[j] = temp * temp;
            }
        }
    }

   public RandomStream getStream () {
        // It is assumed that stream is always the same
        // as the stream in the InverseGaussianGen.
        return stream;
    }

   public void setStream (RandomStream stream) {
       this.stream = stream;
    }\end{hide}

   public int getNumberOfRandomStreams() \begin{hide} {
        return numberOfRandomStreams;
    }\end{hide}
\end{code}
\begin{tabb} Returns the number of random streams of this process.
It is useful because some subclasses use different number
of streams.  It returns 1 for \class{InverseGaussianProcess}.
\end{tabb}
\begin{code}
\begin{hide}
}
\end{hide}
\end{code}
